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Physics has played a fundamental role in medicine sciences, specially in imaging diagnostic. Cur- 
rently, image reconstruction techniques are already taught in Physics courses and there is a growing 
interest in new potential applications. The aim of this paper is to introduce to students the electrical 
impedance tomography, a promising technique in medical imaging. We consider a numerical exam- 
ple which consists in finding the position and size of a non-conductive region inside a conductive 
wire. We review the electric impedance tomography inverse problem modeled by the minimization 
of an error functional. To solve the boundary value problem that arises in the direct problem, we 
use the boundary element method. The simulated annealing algorithm is chosen as the optimization 
method. Numerical tests show the technique is accurate to retrieve the non-conductive inclusion. 
Keywords: electrical impedance tomography; boundary element method; simulated annealing al- 
gorithm; inverse problem; optimization. 

A Hsica tem tido um papel fundamental nas ciencias medicas, especialmente em diagnosticos por 
imagem. Atualmente, as tecnicas de reconstrugao de imagem ja sao ensinadas nos curso de Fi'sica 
e existe um crescente interesse em possfveis novas aplicagoes. O objetivo deste trabalho e apresen- 
tar aos alunos a tomografia de impedancia eletrica, uma promissora tecnica de imageamento em 
medicina. Para isso, consideramos um exemplo numerico que consiste em encontrar a posicao e o 
tamanho de uma regiao nao condutora no interior de um fio condutor. Nos revisamos o problema 
inverso da tomografia de impedancia eletrica modelado pela minimizacao de um funcional de erro. 
Para resolver o problema de valor de contorno que surge no problema direto, nos usamos o metodo 
dos elementos de contorno. O algoritmo de recozimento simulado foi escolhido como metodo de 
otimizacao. Testes numericos mostram que a tecnica e precisa para encontrar a inclusao nao con- 
dutora. 

Palavras-chave: tomografia de impedancia eletrica; metodo dos elementos de contorno; algoritmo 
de recozimento simulado; problema inverso; otimizacao. 



1. INTRODUCTION 

Electrical impedance tomography (EIT) is a technique 
to obtain the inner image of an object exploring its elec- 
trical properties. It consists to apply an electric volt- 
age (or electric current) stimulus through electrodes po- 
sitioned on the object external surface. The correspond- 
ing response (electric current or voltage) is measured by 
the same electrodes. The obtained data are supplied to 
a computer, which has an algorithm to reconstruct the 
conductivity distribution inside the object. This distri- 
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bution can be interpreted as the interior object image. A 
scheme of a medical EIT is depicted in Fig. [TJ 

The first occurrence of imaging the human body in 
terms of its electricalproperties occurred in 1978 by Hen- 
derson and Webster [l| . After this study, several research 
groups have emerged providing significant developments 
in EIT technique. Current applications of EIT in medical 
area are the monitoring of pulmonary functions, breast 
cancer detection, blood flow monitoring inside the heart 
and imaging of human brain functions. Also, in engi- 
neering, EIT applications includes the monitoring of two- 
phase flow in a pipe, detection of corrosion or cracks in- 
side metallic objects and retrieval a pipeline route in the 
subsoil 0-Q 

Although the EIT technique presents attractive fea- 
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FIG. 1. Layout of the electrical impedance tomography. An 
electrodes belt is attached on the body surface. A voltage 
source applies a stimulus on the electrodes and the corre- 
sponding response is measured and supplied to a computer. 
A specific software reconstructs the interior body image. 

tures, such as low cost, portability and robustness, it is 
not widely used yet due to the difficulty to obtain a good 
image resolution. This difficulty arises because the EIT 
image reconstruction is an inverse-problem and therefore, 
it is intrinsically ill-posed Q. The Hadamard criteria, (i) 
existence, (ii) uniqueness and (iii) continuous data depen- 
dence on the solution are not simultaneously guaranteed. 
The first criterion is satisfied, since the physical body be- 
ing imaged certainly has a actual conductivity distribu- 
tion. With respect to the second one, it has been shown 
for the Dirichlet-to-Ncumann map there is a unique solu- 
tion for the conductivity distribution inside the domain 
@. However, the third criterion has been the hardest 
one to be overcome. To deal with it, the researchers have 
tried to find a good and suitable regularization function 
or alternative approaches. For example, it is possible to 
restrain the search-space to speed up the convergence to 
the actual solution and avoid nonphysical solutions Q • 

A typical way to solve the reconstruction EIT problem 
is to treat it as an optimization problem. In this case, one 
must minimize an error functional which is a value ex- 
pressing the discrepancy between two different models of 
the same problem, one experimentally obtained and the 
other one numerically calculated. The conductivity dis- 
tribution that yields the error function global minimum 
corresponds to the sought image. 

Due to the ill-conditioned nature of the inverse EIT 
problem, the optimization surface presents irregularities 
such as multiple local minima and almost flat regions. 
Hence, one must adopt an efficient optimization proce- 
dure to handle such topological features and we have 
chosen the Simulated Annealing (SA) d, Q . This is 
a stochastic method in which the main feature is the 
possibility of eventually accept a solution that increases 
the objective function. This allows the method to detrap 
from local minima basins or flat regions and reach the 
global minimum E| ■ 

To carry out the iterative searching, one must solve the 
direct problem with different conductivity distributions. 
It corresponds to a boundary value problem (BVP), 
which models the stimulus/response process. Frequently, 
this BVP does not have analytical solution, requiring 
the use of numerical methods, such as the Finite Ele- 



ment Method (FEM) and the Finite Difference Method 
(FDM). Here, the Boundary Element Method (BEM) 
was chosen. This is well suited to EIT problem, since 
it requires only the discretization of the domain bound- 
ary and the solution on the boundary is calculated firstly 
without calculating the values of the electric potential 
inside the domain [l2|, [l3j]. Moreover, BEM offers im- 
portant advantages such as great flexibility for arbitrary 
geometries and boundary shape and easiness of imple- 
mentation [3 EH ■ 

Our goal in this paper is to present to students the EIT 
technique using BEM and SA algorithm. This study was 
motivated by the difficulty in teaching the physical and 
mathematical methods/processes related to image recon- 
struction. Some studies have been carried in experimen- 
tal approach, such as the assembly proposed by Mylott et 
al. [Ty] to study the computed tomography. Considering 
the students' growing interest in computational physics 
in the last years, we have proposed a numerical example 
to introduce the EIT technique. It consists in retrieving 
a non-conductive region inside a conductive wire. In en- 
gineering, this non-conductive region may represent, for 
example, a manufacturing defect or a normal wear in the 
wire which could compromise its use in an apparatus. In 
medicine, it may model an air bubble inside an artery 
hampering the blood flow. 

Our numerical example considers a conductive cylin- 
drical wire of length 10.0 and diameter 1.0 (arbitrary 
units). The non-conductive region is a sphere of radius 
0.3 with center located along the cylinder axis at 7.0 units 
far from one of its extremities, as shown in Fig. [2j The 
challenge is summarized as a two variables optimization 
problem, the center ^-coordinate and radius R of the 
non-conductive region. Due to the domain symmetry, 
we have taken a longitudinal section of the wire as a 
two-dimensional domain and placed the left bottom cor- 
ner at the origin of a Cartesian system xOy. Also, the 
stimulus/response process consists in applying an elec- 
tric voltage V a b = 12 V between the wire extremities a 
and b (V a = 12 V and V& = V) and measuring the cur- 
rent flux on the same extremities. Moreover, the lateral 
area of the wire is considered electrically insulated, such 
that the current flux vanishes. 

insulating 







7.0 
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FIG. 2. Scheme of the numerical example: cylindrical con- 
ductive wire with length 10.0 and diameter 1.0 with a non- 
conductive region of radius 0.3 located 7.0 units far from ex- 
tremity a. 

In Sec. [3J we review the EIT mathematical modeling 
based on the functional approach, the Boundary Ele- 
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merit Method and the Simulated Annealing algorithm. 
In Sec. [3J we show the numerical tests and the results. 
Finally, in Sec. IU we discuss the results and conclude. 

2. PROBLEM STATEMENT AND METHODS 

Here, we present the mathematical and numerical ba- 
sis of EIT problem, specially focused on our numerical 
example. 

I. EIT mathematical formulation 

The mathematical modeling of the EIT problem is ob- 
tained from the Maxwell's equations. Since the stimulus 
processing is made at low frequency, the inductive ef- 
fects can be ignored. Also, considering the domain fully 
ohmic, we can neglect the capacitive effects. The elec- 
tric field is E = — V0, where <fi is the electrical potential 
and V is the gradient operator. The current density is 
J = —oE = —aS7(j), being a the electrical conductivity. 
Considering no internal current sources, V • J = and 
the partial differential equation (PDE), 

V • (-aV<t>} = in ft, (1) 

governs the electrical potential <p inside the domain ft. 
To simulate the stimulus/response process, we have two 

boundary conditions given by <p = 4> an( l — °~~zr = J ■> 

on 

where <f> ls the voltage profile and J is the normal com- 
ponent of the electric current flux, both on the boundary 
dft. 

The numerical example considers a two-phase medium 
problem. One of them (the spherical non-conductive re- 
gion) has null-conductivity and the other one (the rest of 
the wire) has constant and uniform non-vanishing con- 
ductivity. In this case, the interior of the non-conductive 
inclusion is interpreted as the exterior of the domain and 
the conductivity distribution a corresponds to the posi- 
tion and size of the non-conductive region. Also, Eq. ([T]) 
is simplified to the Laplace equation 



v2 ^ = ^ + ^ = < for (*>»)en. (2) 

Considering unitary the conductivity of the non null 
phases (a = 1), the boundary conditions become 



The stimulus-response can be mathematically mod- 
eled through the Dirichlet-to-Neumann map, in which a 
known voltage profile 4> is imposed on the boundary and 
the electric current fluxes J are calculated on the same 
boundary. Otherwise, if a known electric current flux 
profile J is imposed and the voltage cj> is calculated, the 
model is called Neumann-to-Dirichlet map. Although in 
practice the EIT problem is formulated as a Dirichlct- 
to-Ncumann map, mathematically, it represents a mixed 
problem since it is imposed null current flux for the in- 
ternal boundary, i.e., Eq. @ vanishes. 

The direct problem consists in solving the BVP given 
by Eqs. © and (J3J, considering a knowing conductiv- 
ity distribution a to obtain the current flux from Eq. Q . 
However, in EIT image reconstruction, a is unknown and 
we have an inverse problem. One approach to solve it 
and to obtain a is to construct and minimize an error 
functional (objective function) that compares the mea- 
surements obtained from two different models, the actual 
model and the numerical model. 

For the Dirichlet-Neumann map, the actual model cor- 
responds to a collection of measurements, 3 actual, which 
are the current fluxes on the electrodes obtained from ex- 
perimental assembly. This model contains the unknown 
actual conductivity distribution, o- ac tuai, which must be 
found. The numerical model considers a known prospec- 
tive conductivity distribution, a prosp , to solve the direct 
problem to obtain 3 num , which corresponds to the nu- 
merical current fluxes on the same electrodes. The error 
functional must be defined as a discrepancy measure be- 
tween 3 actual and 3 num . Here, we have defined the error 
functional e{a prosp ) as the mean square error between 
the experimental and numerical current fluxes, 



e(o~prosp) — / j \ J actual ^num\ ' (5) 



where J„ ctua ; and Jnum are, respectively, the actual and 
numerical fluxes in the electrode i and n is the number 
of electrodes. 

The idea of the optimization process is to perform an 
iterative search to find the prospective conductivity dis- 
tribution that minimizes Eq. ([5]). In each iteration, a 
prospective solution, a prosp , is generated and supplied to 
Eqs. (|2|) and ([3]) to yield a numerical flux 3 num . This flux 
is compared with the actual one 3 actual to calculate the 
error functional, from Eq. ([5]). This process is repeated 
so that e {o- prosp ) « and, consequently, a prosp m a actua i 
within an acceptable error level. 



<j)(x, y) = <j)(x, y), for (x, y) E d£l, (3) II. Boundary element method 

The essence of BEM is to convert the field equation to 
d(f>(x, y) — integral equations on the domain boundary through the 

qZ = J( x >y)> f° r [x,y) € dO,. (4) reciprocal relation, 
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(6) 



where <pi and 02 are two general solutions of Eq. ^ , and 
through the fundamental solution, 



*(x t y; i,r{) = — ln[(z - £) 2 + (2/ - '/) 2 L (7) 

which is a particular solution of Eq. $Z§. From Eqs. © 
and (O, it is possible to obtain two integral equations, 
one for the points (£, 77) inside the domain and an- 
other for the points (£,77) on the boundary dCl. In the 
discrctized form, these integral equations are 



fe=i 



(8) 



for (£, r?) £ f2, and 



1^ 

2 



fori) = E [W«.1)-^ W «.1)] - (9) 



k=l 



for (£,77) G dfl, where 



(10) 



— [$(x, 7/; £,77)] ds(x,y) (11) 



The external boundary is discretized in iV ex i straight 
elements, dVt^ k \k = 1,2, N ext , and the internal one 
in N lnt elements, dfl^,k = N ext + l,N ext + 2, ...,N T . 
The coordinates (x^ k ',y^) of the elements extremities 
are supplied to the program and, by convention, the el- 
ements are numbered following the counter clockwise di- 
rection for the external boundary and clockwise for the 
internal one. This guarantees the normal unitary vector 
n always points to the outside of the domain. Also, for 
each boundary (external and internal), the final extrem- 
ity of the last element is made to coincide with the initial 
extremity of the first one. In this case, the domain and 
the inclusion are approximated by polygons. 

Equation @ produces a linear system of Nt equations 
with Nt unknowns, either <p or J. Solving this system 
provides the pair (0, J) for all boundary elements. Hence, 
the direct problem solution of EIT is completed. If the 
solution <j) at some internal point is required, it can be 
calculated through Eq. ©. 



III. Simulated annealing algorithm 

The simulated annealing algorithm is a stochastic op- 
timization method that has its origins in the Metropolis 
acceptance criterion when two system configurations are 
compared [l7|. It is based on the probability to find the 
system with energy E, given by the Boltzmann weight, 
e -E/k B T ^ wnere j> is the temperature and ks is the Boltz- 
mann constant. 

In 1983, Kirkpatrick et al. [H| solved the salesman 
problem using the Metropolis criterion adding an impor- 
tant differential. They adopted a cooling schedule for 
the temperature to control the search stochasticity. In 
the thermodynamic framework, if the temperature of a 
liquid material is slowly cooled down, the atoms arrange 
themselves to form a structure (a perfect crystal) , which 
has the lowest internal energy state. However, if the cool- 
ing process is not sufficiently slow, the final structure is 
not perfect and the internal energy is not the lowest one. 

The analogy with an optimization process is made con- 
sidering the objective function f(x) as the energy E of 
the system, the solutions the system configurations 
and the temperature T becomes a control parameter of 
the process. The optimization process is done iteratively 
starting from an initial solution xo and temperature To- 
At each iteration, k searches are performed. Each new 
solution, Xi+\ (i = 0, 1, 2, k), is generated through a 
pre-defined visitation distribution and compared with the 
current one, x,, to be, or not, accepted according to the 
Metropolis criterion. If /(xj+i) < f(xi), the new solu- 
tion is accepted and if /(xi+i) > f(xi), the new solution 
is accepted with probability 



(12) 



where A/ = f(x i+ i) — f{x*i) and ks = 1- After the 
last solution Xk is evaluated, the temperature is de- 
creased through a cooling schedule and the search process 
restarts. The last solution accepted at the previous tem- 
perature is taken as the initial solution to the current 
one. 

The most common ways to generate a new solution 
Xi + i from the previous one Xi are through uniform or 
Gaussian deviates. Here, we have chosen the Gaussian 
one in the following form 



To 



A 77, 



(13) 



where To is the initial temperature, T t is the tempera- 
ture at the iteration t, A is a parameter that depends on 
the search space range and 77 is a random number with 
normal distribution N(0, 1). 

According to Eqs. (fl~2j) and ([T3]). at high temperatures, 
solutions can be generated relatively far from the previ- 
ous one, characterizing a global search. Also, there is 
high probability to accept a new solution that increases 
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the objective function. As the temperature decreases, the 
search becomes essentially local and the acceptance of an 
uphill error solution becomes less likely to occur. Hence, 
the temperature T plays a fundamental role. The choice 
of a suitable cooling schedule is essential to ensure a good 
performance of the optimization process, i.e., to reach the 
global minimum of the objective function with as few 
steps as possible. The most commons cooling schedules 
in the literature are the logarithmic and the geometric 
ones. We have chosen the last one in which the temper- 
ature T decreases with the iteration t according to 

T t = a t T , (14) 

where To is the initial temperature and a (E (0, 1) is the 
cooling rate. 

3. NUMERICAL TESTS AND RESULTS 

To solve the numerical example proposed in Sec. [TJ 
we developed in MATLAB© lan guage a program that 
includes the BEM and the SA algorithm routines. 

To construct the error functional, Eq. one must 
confront two models of the same problem, the actual and 
the numerical ones. The actual model is simulated nu- 
merically, since it is not the goal of this paper to explore 
the experimental part of the problem. Hence, the ac- 
tual current flux, 3 actual , is obtained solving numerically 
the Eq. (J5|) through BEM with the boundary conditions, 
0(0, y) = 12V and (10,t/) = V, for y e [0,1], and 
J(x,0) = J(x,l) = 0, for x e [0,10]. Also, it was 
adopted J (x, y) = for the points (x, y) on the bound- 
ary of the non-conductive spherical region. The numeri- 
cal model is defined as the actual one. However, its non- 
conductive region, called prospective inclusion, has center 
x-coordinate (x p ) and radius (R p ) variables in order to 
carry out the searching process. 

For both models, the external boundary discretization 
was made using 220 elements, being 100 elements for each 
side, top and bottom, and 10 elements for each side, left 
and right. The internal inclusion was discretized in 80 
elements, totalizing 300 boundary elements. 

I. Test 1 - Error functional behavior 

The first test explores the behavior of the error func- 
tional surface. The variables x p and R p were systemati- 
cally modified over the search space [2.0, 8.0] x [0.1, 0.4] in 
a regular mesh of 61 x 61 = 3721 points. The error func- 
tional is calculated for each position x p and radius R p 
generating the error functional surface, as shown in the 
plot of Fig. [3] with error functional axis in logarithm scale. 
It is possible to see almost flat regions and a narrow chan- 
nel that contains a prominent global minimum. These 
characteristics require a powerful optimization method 
to reach the global minimum. 



error functional surface 
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FIG. 3. Error functional surface in logarithm scale. It 
presents almost flat regions and a channel that contains a 
prominent global minimum. 



II. Test 2 - Applying the Simulated annealing 

In this test, we evaluated the performance of the devel- 
oped program to retrieve the position and radius of the 
non-conductive region. Since SA algorithm is stochastic, 
its convergence to the actual solution is not always guar- 
anteed. Therefore, it is necessary to run the program 
many times to assess its performance. Also, to make 
a more careful analysis, we have divided the Test 2 in 
three parts. 

In the first part, Test 2. a, the intent was to find only 
the center ^-coordinate of the non-conductive region. We 
set the prospective inclusion radius equal to the actual 
one (R p = 0.3). The initial solution of x p was generated 
randomly in the range [1.0,9.0] and, to generate a new 
solution x p k+1 ^ from the previous one x p by Eq. (fl"3j) . 
it was chosen A = 0.8. In the second one, Test 2.b, the 
program searched only the non-conductive region radius. 
In this case, we set the prospective inclusion center x- 
coordinate equal to the actual one (x p = 7.0). The initial 
solution of R p was generated randomly, but in the range 
[0.05,0.45]. The new solution R p k+1 ^ from the previous 
one R p k ^ by Eq. ifT3")) was generated considering A = 0.04. 
For both, Test 2. a and Test 2.b, the program was exe- 
cuted considering 1000 iterations, with k = 1, and we set 
a = 0.95 and T a = 1000. 

Finally, the third part, Test 2.c, considers the search 
for both, position and radius. In this case, we joined the 
two programs used in the previous parts in a single one. 
More specifically, in each iteration (at each temperature) 
the program carried out two separated searching, one 
for the position and another for the radius. Due the 
difficulty to carry out the optimization searching with 
two variables, we adopted 2000 iterations, with k = 1, 
and a = 0.97. As the previous parts, we set To = 1000, 
A = 0.8, to generate a new x p , and A = 0.04, to generate 
a new R p . Again, the initial solutions of x p and R p were 
generated randomly in the same ranges of the previous 
parts. 
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propective inclusion radius X iteration 




For each part of the Test 2, the program was exe- 
cuted 50 times and the final solution x p and/or R p found 
in each one was stored. Table U shows the minimum, 
maximum, mean and standard deviation of the final so- 
lutions considering all of the runnings for each part. For 
the best running of the Test 2.c, the solution found was 
Xp = 7.0005 and R p = 0.3000, with error functional of 
2.6 x 10~ 17 . For the best running, the evolution of the 
error functional, the prospective center x-coordinate and 
the prospective radius during the optimization process 
are shown in Fig. 



TABLE I. Minimum, maximum, mean and standard deviation 
of the final solutions x p and/or R p for the three parts of the 
Test 2, considering 50 runnings. 





Test 2. a 


Test 2.b 


Test 2.c 






R P 




R p 


min 


6.9925 


0.2995 


6.8388 


0.2445 


max 


7.0081 


0.3005 


7.5467 


0.3163 


mean 


7.0000 


0.3001 


7.1203 


0.2876 


std 


0.0032 


0.0002 


0.2053 


0.0214 



(c) 



4. DISCUSSION AND CONCLUSION 

We have presented to students the electrical impedance 
tomography technique and the application of retrieving 
the size and position of an non-conductive inclusion in- 
side of a conductive wire using the boundary element 
method and simulated annealing algorithm. Numerical 
simulations assessed either the behavior of the error func- 
tional in the search-space and the performance of the 
implemented program. The results show that the prob- 
lem has a difficult error surface to be optimized. Despite 
this, the implemented program has presented a satisfac- 
tory accuracy to find the non-conductive region. Finally, 
the presented algorithms and methods have a wide scope 
of application and should be better appreciated in topics 
of undergraduate and graduate physics programs. 



ACKNOWLEDGMENTS 



FIG. 4. Evolution with the iteration for the best running of 
Test 2.c of the (a) error functional (b) prospective center 
i-coordinate and (c) prospective radius. 



ASM acknowledges Conselho Nacional de Desenvolvi- 
mento Cienti'fico e Tecnologico (CNPq) (305738/2010-0). 
VR acknowledges Fundagao de Amparo a Pesquisa do 
Estado de Sao Paulo (FAPESP) (2008/01284-4). 



[1] R. P. Henderson and J. G. Webster, IEEE Transactions 
on Biomedical Engineering 25, 250 (1978) 

[2] W. Lionheart, N. Polydorides, and A. Borsic, "Electrical 
impedance tomography: methods, history, and applica- 
tions," (Institute of Physics Publishing, 2005) Chap. The 



reconstruction problem 
[3] J. F. Abascal, S. R. Arridge, D. Atkinson, R. Horesh, 
L. Fabrizi, M. D. Luca, L. Horesh, R. H. Bayford, and 
D. S. Holder, Neurolmage 43, 258 (2008) 



7 



[4] A. Osella, G. Chao, and F. Sanchez, American Journal 
of Physics 69, 455 (2001) 

[5] L. Borcea, Inverse Problems 18, R99 (2002) 

[6] A. I. Nachman, The Annals of Mathematics, Second Se- 
ries 143, 71 (1996) 

[7] O. Menin and V. Rolnik, International Journal of Modern 
Physics C 22, 825 (2011) 

[8] E. Aarts and J. Korst, Simulated Annealing and Boltz- 
mann Machines: A Stochastic Approach to Combinato- 
rial Optimization and Neural Computing (J. Wiley & 
Sons, 1988) 

[9] D. Goldberg, Genetic Algorithms in Search, Optimiza- 
tion and Machine Learning (Addison- Wesley Longman 
Publish, Boston, 1989) 

[10] T. C. Martins, E. D. L. B. Camargo, R. G. Lima, M. B. P. 
Amato, and M. S. G. Tsuzuki, in Engineering in Medicine 
and Biology Society, EMBC, 2011 Annual International 
Conference of the IEEE (2011) 

[11] H. Kim, C. Boo, and Y. Lee, International Journal of 
Control, Automation and System 2, 211 (2005) 



[12] V. Rolnik, O. H. Menin, G. L. C. Carosio, and P. J. 
Seleghim, in ABCM Symposium Series in Mechatronics, 
Vol. 4 (2010) pp. 835-843 

[13] R. Duraiswami, K. Sarkar, and G. L. Chanine, Engineer- 
ing Analysis with Boundary Elements, 13(1998) 

[14] C. A. Brebbia and J. Dominguez, Boundary Elements 
Methods An Introductory Course, 2nd ed. (Compu- 
tational Mechanics Publications Southampton Boston, 
1992) 

[15] W. T. Ang, A Beginner's Course in Boundary Element 
Methods (Universal Publishers, Boca Raton, Florida, 
2007) 

[16] E. Mylott, R. Klepetka, J. C. Dunlap, and R. Widenhorn, 
European Journal of Physics 32, 1227 (2011) 

[17] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, 
and E. Teller, Journal of Chemistry Physics 21, 1087 
(1953) 

[18] S. Kirkpatrick, C. Gellat, and M. Vecchi, Science 220, 
671 (1983) 



